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A  new  iterative  method  is  presented  for  solvinq  non-symmetric  linear 
systems  of  equations.  The  method  requires  that  the  symmetric  part  of  the 
matrix  of  the  linear  system  be  positive  definite,  and  the  method  is  efficient 
only  if  the  symmetric  part  is  easily  invertible.  The  method  is  modeled  on  the 
conjuqate  qradient  method  for  symmetric  positive  definite  systems  and  has  the 
finite  termination  property.  The  results  from  several  numerical  experiments 
are  presented  and  compared  with  a  similar  method  proposed  by  Concus,  Golub, 
and  Widlund. 
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SIGNIFICANCE  AND  EXPLANATION 


In  this  report  a  new  method  is  presented  for  the  solution  of  linear 
systems  of  equations. 

Ax  =  b  , 

where  the  matrix  A  is  a  non-symmetric  matrix.  The  matrix  A  can  be  written 
as  the  sum  of  its  symmetric  and  skew-symmetric  parts 

A  *  P  +  Q 

and  the  method  requires  that  the  symmetric  part ,  P,  be  positive  definite, 
i .  e . 

( x , Px )  >  0 

for  all  non-zero  vectors  x  .  the  method  is  efficient  only  when  the  solution 
of  the  linear  system  Py  ■  c  can  be  obtained  easily.  This  is  the  case  in 
many  problems  such  as  the  solution  of  elliptic  equations  whose  hiqhest  order 
part  is  the  Laplacian. 

The  method  is  modeled  on  the  conjuqate  gradient  method  which  is  a  widely 


used  method  for  symmetric  positive  definite  systems. 


lhe  responsibility  for  the  wordinq  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


A  GENERALIZED  CONJUGATE  GRADIENT  METHOD 
FOR  NON-SYMMETRIC  SYSTEMS  OF  LINEAR  EQUATIONS 

John  C.  Rtrikwerda 

1.  Introduction 

In  this  paper  a  generalized  conjugate  gradient  method  for  solving  linear 
systems  of  equations  is  presented.  The  proto-type  for  this  method  is  the 
system 

(1.1)  (I  +  S)x  =  h 

where  S  is  a  skew-symmetric  matrix  and  I  is  the  identity  matrix.  The 
method  is  derived  as  an  acceleration  of  the  steepest  descent  method  for  the 
system  (1.1)  where  the  norm  of  the  residual  is  minimized. 

Concus  and  Golub  (1976)  and  Widlund  (1978)  have  presented  and  discussed  a 
generalized  conjugate  qradient  method  for  non-symmetric  linear  systems  which 
has  some  similarities  with  the  present  method  (see  also  Hageman  et  al. 

1980).  Althouqh  hoth  methods  can  he  derived  in  several  ways,  the  derivation 
of  the  method  of  this  paper  is  different  in  spirit  from  the  method  of  Concus 
and  Golub  (1976)  and  Widlund  (1978)  as  qiven  in  their  papers.  Their  method  is 
derived  by  imposing  orthogonality  constraints  on  an  appropriate  sequence  of 
vectors  and  the  convergence  properties  are  then  deduced.  The  method  presented 
here  is  derived  as  an  acceleration  of  a  steepest  descent  method  and  the 
orthoqonality  results  then  follow.  As  shown  in  section  A  the  two  methods  have 
similar  rates  of  convergence  and  in  the  numerical  experiments  they  behaved 
similarly. 

The  method  of  this  paper  can  he  applied  in  Hilbert  spaces  as  was  done  by 
Widlund  (1978),  hut  little  would  be  gained  by  the  extra  generality  so  we  will 
consider  the  method  onlv  for  finite-dimensional  spaces. 
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2.  Derivation  of  the  Method 


Consider  the  system  of  linear  equations 

(2.1)  Ax  -  b 

where  A  is  a  real  n  *  n  matrix  and  b  is  a  real  n-vector.  We  decompose 
A  into  its  symmetric  and  skew-symmetric  parts, 

(2.2)  A  *=  P  +  p 

where  p  is  the  symmetric  part  of  A  and  P  is  the  skew-symmetric  part, 
i.e. 

1  T  1  T 

P=j(A+A  ),  P  =  —  (A  -  A  )  . 

We  assume  that  P  is  positive  definite  and  hence  the  system  (2.1)  has  a 
unique  solution. 

We  will  begin  by  considering  the  special  case  of  equation  (1.1)  where 
P  is  the  identity  matrix.  In  section  5  we  will  show  that  the  method  can 
treat  those  cases  where  the  system 

Pz  =  c 

can  be  easily  solved.  This  is  the  same  situation  considered  by  Concus  and 
Golub  (1976)  and  Widlund  (1978). 

The  standard  conjugate  gradient  method  (Hestenes  and  Stiefel  (1952))  is 
used  to  solve  linear  systems  such  as  (2.1)  when  the  matrix  A  is  symmetric 
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(2.4) 


Yk  =  ^  /(pk»APk) 


5k  =  |rk+1|2/|rk|2  • 


By  analogy  with  the  conjuqate  qradient  method  for  positive  definite 

systems  we  consider  the  following  iterative  scheme  for  the  system  (1.1) 

,  k+1  k  k 

a)  w  =  w  +  a^p 


(2.5) 


.  .  k+1  k  k 

b)  r  =  r  a^d+SJp 

.  k  +  1  k+1  -  k 

c)  p  =  r  -  @kp 


The  parameters  «k  and  8k  are  to  be  chosen  to  minimize  |rlc+1|2 
k  k  1 

given  r  and  p  .  Equation  (2.5b)  is  a  consequence  of  the  definition  of 
the  residual  vector,  rk  =  b  -  (I+Sjw*,  and  (2.5a). 

One  obtains  a  steepest  descent  method  by  settinq  all  3  =  0,  so 

V  V 

p  =  r  ,  and  then  ak  =  a^,  where 

-  |rk|2/(|rk|2  +  |spk|2) 

minimizes  |rk+1|  . 

i  k + 1  i  ^  k 

To  derive  our  method  we  choose  ak  to  minimize  |r  |  given  p  .  We 


(2.6) 


|rk+1|  =  |rk|  -  »k(rk,  (I+S)pk)  +  c*k| (I+S)pk| 2 


and  so 


(2.7) 


<*k  =  (rk,  (I+S)pk)/(  | pk |  2  +  |spk|2) 


is  the  optimal  value  of  ak>  Ttie  first  consequence  of  (2.7)  is  that 


(2.8) 


(rk  +  \  ( I+S )pk )  =  0 


by  (2.5b),  and  secondly,  usinq  (2.8)  with  (2.5c) 

(rk,  ( I+S)pk)  =  (r\  (I+S)(rk  -3k_1pk"1)) 

=  (rk,  (l+S)rk)  -3.  ,(rk,  (I+S)pk_1) 

2 

=  Irk!  . 


Hence 


(2.9) 


ak  *  |rk|  /( |pK| 2  +  |Spk|  2)  . 
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The  relation  (2.6)  then  becomes 


(2.10)  )r  | 

and  since  |rk+1 |  <  |rk|  we  have 


=  r  (1  -  a  ) 


0  <  <X  <  1  . 


From  (2.R)  and  (2.5b)  it  also  follows  that 

2 

, ,  .  i  k+1 1  .  k+1 

(2.11)  I r  I  *  (r 


,  k+1  k. 

(r  ,  r  )  . 


We  now  determine  8^  ^  piven  r  and  p  .  We  see  from  (2.10)  and 


k  i_k+1 


(2.9)  that,  piven  r  ,  |r 


is  minimized  when  a,  is  maximized  and  this 


repuires  that  |p*|  +  |Sp  |  he  minimized.  We  have 


|Pk|  +  |spk|  =  ( ( I+S )pk,  (I+S)pk) 


-  | (I+S)rk |  -  28k_1( (I+S)rk,  (I+S)pk_1) 


+  ^,l(i+s)pk'1' 


and  hence 


(2.12) 


2  2 

8k_1  *  (d+S)rk,  (I+S)pk_1  )/( |pk_1  |  +  |spk_1  f  )  . 


From  (2.12)  and  (2.5c)  we  have 


(2.13) 


((I+S)pk,  (I+S)pk_1)  =  0  , 


and  hv  subtractinp  times  (2.13)  from  (2.R)  with  k  -  1  replacinp  k 

we  obtain 


(2.14) 


,  k+1  ,,  k-1 

(r  ,  ( I+S )p  )  =  0 


We  now  wish  to  obtain  an  expression  for  8^  which  is  simpler  than 


(2.12).  We  have 


|f  ±1  V  Ifxl  V 

((r+S)r  ,  (I+S)p  )  =  (Sr  ,  (I+S)p  )  hv  (2.R) 


1  k+1  k  k  +  1 .  ,  _  . 

=  —  (Sr  ,  r  -  r  )  by  (2.5b) 

1  k+1  k 

■%,sr  'r’  ■ 


Hence 
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2 


a  k+1  k .  |  k  | 

8k  =  (Sr  ,  r  )/ | r  | 


However, 


k+1  k  k+1  _  k. 

(Sr  ,  r  )  =  -(r  ,  Sr  ) 


-(rk+1,  (I+S)rk)  +  |r)C+1|  by  (2.11) 


,  *+1  k  o  k-1 . .  |  k+1  | 

-(r  ,  (I+S)(p  +  3k_1p  ))  +  |r  |  by  ( 


1c  ^  1 

=  |r  |  by  (2.8)  and  (2.13)  . 


So 

(2.15) 


2  2 

0k  =■  |r  |  /|  r  |  -  1  -  ak  by  (2.10)  . 


We  now  summarize  the  alqorithm. 
(2.16) 


.  k+1  k  k 

a)  w  *  w  +  <*kp 


wi  th 


and  8,  =  1  -  a,  . 
k  k 


,  .  k+1  k  _  k 

b)  r  ■  r  -  ak(I+S)p 

,  k+1  k+1  0  k 

c)  p  -  r  -  8kp 


d)  ay  -  |rk|  /(|pk|  +  |Spk|  ) 


For  initial  vectors  we  take  w°  arbitrary,  ru  =  p 


.5c) 


I 


-  ( I+S )w°. 


3.  Orthoqonalitv  Relations. 

The  purpose  of  this  section  is  to  prove  orthoqonalitv  relations  for  the 
vectors  qenerated  by  the  ahove  alqorithm.  The  main  result  is  contained  in  the 
followinq  theorem. 

Theorem  3. 1 

For  the  alqorithm  (2.16) 

(3.1)  ((I+SJp1,  (I+S)pj)  =  0 
for  i  jt  j . 

Proof 

We  beqin  by  obtaininq  a  three  term  recurrence  relation  for  the  vectors 
k  k 

p  .  By  eliminatinq  the  vectors  r  from  (2.16b)  and  (2.16c)  we  obtain 

. ,  k+1  _  k  k-1 

(3.2)  p  +  a  Sp  -  =  0 

for  k  >  0  where  p-1  *  0. 

First  we  show  that  (3.1)  is  true  for  i  =  1 ,  i  =  0.  This  is  immediate 
from  (3.2)  and  the  skew-svmmetrv  of  S 

((I+SJp1,  (I+S)p°)  =  -ao((I+S)SpP,  (I+S)pn) 

=  0  . 

Now  assume  that  (3.1)  is  true  for  k  >  i  >  i  >  0,  we  will  show  that  it  is 
also  true  for  k+1=i>j>0.  If  i=  k+1  and  j  =  k,  (3.1)  is  true  by 
(2.13).  If  i  =  k+1  and  j  =  k-1,  we  have  by  (3.2) 

((I+S)pk+1,  (I+S)pk-1) 

(3.3)  .2 
=  -o^  ( (I+S)Sp  ,  (I+Sip"1-')  +  Pk_1|(I+S)p  '  | 

Consider  the  first  expression  on  the  riqht-hand  side  of  (3.3). 

<{I+S)Spk,  (I+S)pk_1) 

-  - ( ( I +s ) pk ,  (I+S)Spk_1) 


-6- 


' 


-  -((I+S)pk,  <i+s)(-pk  +  Pk"2))/ak"1 

2 

-  l(I+S)p*|  /c«k_i  . 

Thus  (3.3)  is  equal  to 

a  2  2 

-- —  |(I+S)pk|  +  0  |(I+S)pk‘1  I  , 

k-1  1 


and  by  the  expressions  (2.16d)  and  (2.15)  this  is  equal  to 


I  (I+S)pk_1  |  +  |  (I+S)Pk_1 


k-1 


|rk-1|' 


-  0  . 

Now  for  i  -  k+1  and  j  <  k-1  the  relation  (3.2)  follows  easily. 
(<X+8)pk+1,  (I+S)pj) 

-  -ak((l+S)Spk,  (I+S)pj)  +  0k_,<  nr+s>Pk*  (I+S)pj) 

-  c*k((I+S)pk,  (I+S)Spi) 

=  ~  (d+s)Pk,  (i+s)(pj  -  Bj_1Pj"1)) 


*  0  • 


Thus  the  theorem  is  proved • 

V 

It  follows  that  the  vectors  p  are  linearly  independent,  as 
they  are  non-zero,  and  thus  the  algorithm  must  converqe  in  at  most 
where  N  is  the  dimension  of  the  vector  space. 

Other  orthoqonality  relations  are  qiven  in  the  next  theorem. 
Theorem  3. 2 

For  the  algorithm  (2.16) 


(3.4) 

(rk,  (I+S)^)  =  0 

for 

i 

<  k 

(3.5) 

(rk,  <I+S)pj)  -  0 

for 

i 

<  k 

(3.6) 

(p  ,  p  )  =  0 

for 

0 

4  k 

(3.7) 

,  k  j,  ,  k  0 

(r  ,  rJ)  *  (r  ,  r  ) 

for 

1 

<  k 

lonq  as 
N  steps. 
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Proof  of  (3.4) 


For  k  =  1  and  j  =  0  (3.4)  follows  from  (2.8)  since  p°  =  r°.  When 

j  =  k-1  >  0,  (rk,  ( I+S)rk-1 )  =  (rk,  (I+S)(pk_1  +  Bk_2pk“2)  =  0  by  (2.8)  and 
(2.14).  For  j  <  k-1  the  result  follows  from  Theorem  3.1  by  induction,  we 
have 

(rk,  (I+S)^)  =  (rk~\  (I+S)rj)  -  <*k_.,(  (I+S)pk"\  (I+S)rj) 

=  -ak_1((I+S)pk“\  (I+S)(p^+  ®j_lP^_1)) 

=  0  . 

Proof  of  (3.5) 


By  repeated  use  of  (2.16c)  and  (3.4) 

(rk,  (I+S)pj)  =  (rk,  (I+S)rj)  -  0  (rk,  (I+S)pj-1) 

*  -B^lr*,  (I+S)pj-1) 

i  k  0 

*  (-r  n  B.  (r  ,  (I+S)PU)  •  0 

1  =  0  1 

by  (3.4)  since  p°  =  r°. 

Proof  of  (3.6). 


By  (3.2) 


,  k+1  k.  n  .  k-1  k. 
(P  *  P  )  =  6k_.,fP  »  P  > 

k-1 

-  n  6.  (p\p°) 
1=0 


Now  by  (2.16c)  and  (2.16b) 


.1  0.  ,1  0.  Q  »  0 | 
(p  ,  p  )  =  (r  ,  P  )  -  B  IP  | 


(r°,  p°)  -a0((I+S)p°,  p°)  -  6 0 1 p° | 


(r°,  p°)  -  (ao-t60)|p°|  =  0  , 


since  p°  =  r°  and  +  B^  =  1. 


Proof  of  (3.7) 


By  (2.16b)  for  j  <  k 

.  k  j ,  ,  k  j-1,  ,  k  . „  _ .  j-1 . 

(r  ,  r  )  =  (r  ,  rJ  )  -  a^fr  *  (I+S)pJ  ) 

=  (rk,  r1"1)  by  (3.5) 
k  0 

=  (r  ,  r  ),  by  repetition. 

This  proves  all  the  assertions  of  Theorem  3.2. 

V 

The  relation  (3.7)  has  the  geometric  interpretation  that  r  is  on  the 
sphere  of  radius  -j  |r^|  centered  at  —  r^  for  each  j  less  than  k.  This 
again  demonstrates  the  finite  termination  property  of  the  method  since  N 
distinct  spheres  through  the  origin  in  N-space  can  have  only  the  origin  as  a 
common  point. 


_g_ 


4.  The  Rate  of  Convergence 


Considering  the  method  (2.16)  as  an  iterative  method  for  solving  (1.1), 
it  is  natural  to  estimate  the  rate  of  convergence  of  the  method  in  terms  of 
the  spectral  radius  of  S.  For  linear  systems  with  a  large  number  of 
unknowns,  such  as  arise  from  numerical  approximations  to  partial  differential 
equations,  the  finite  termination  property  is  of  little  interest  compared  with 
the  rate  of  converqence. 

Our  first  convergence  results  are  stated  in  the  following  theorem. 

Theorem  4.1 


(4.1) 

(4.2) 


For  the  method  (2.16)  the  following  estimates  hold. 

|rk*,|/|rl'l  <  A 


|rl”2|/|rlt|  <  ( 


/T777 


where  A  is  the  spectral  radius  of  S, 
Proof 


k1 


We  begin  with 

k+1 | 2  ,  k+1  k 

I  =  (r  ,  r  ) 

,  k+1  _  k, 

=  -(r  ,  Sr  ) 

=  ak((l+S)pk,  Srk ) 


by  (3.7) 
by  (3.4) 
by  (2.16b) 


o^( (I+S)rk,  Srk)  -  ak6k_i< (I+S)pk_1,  Srk)  by  (2.16c) 


k,2  ak8k-1  ,  k-l  k  „_k, 
«klSr  I  “  “ -  (r  -  r  ,  Sr  ) 


ak  Sr 


k  .2 


k-1 
a  6 

Itk, 


by  (2.16b) 


by  (3.4)  and  (3.7) 


k-1 


This  gives  the  estimate 

(4.3) 


°k-1  °k 


Sr 


k,  2 


<  A  ‘ 
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since  8k  -  |rk+1 |  /|rk|  . 

The  estimate  (4.1)  follows  from  (4.3),  since 

\  2 
<  A  2 

°k 

and  a  -  1  -  8,  .  So 
k  k 

|rk+1|2/|rk|2  «  Bv  <  A 


k  .2 

1  +  A 


which  is  (4.1). 

The  estimate  (4.2)  is  obtained  by  finding  the  maximum  of  the  product 
ek+16k  *  (1  "°k+1,(1  “  “k*  sub3ect  only  to 


(4.4) 


1 _  1  ,  .2 

- -  +  — —  <  A  +2 

CX  (X 

k+1  k 


which  is  equivalent  to  (4.3).  The  maximum  of  (1  -  ak  +  1Hl  -  a^)  obviously 
occurs  when  equality  holds  in  (4.4).  Thus 

<’  -S.i"1  -  V  -  '  -  Vi  ♦VA 
■  ’  -  ‘V.Vr-  *  3T7  -  ’> 

k  k-1 

*  1  -°k+iVA2  +  15 

and  this  quantity  is  maximized  subject  to  (4.4)  when 

“kli  'T  (,J  ♦  2K 

Bk*i9k‘  ’  -  «2  -  ”/<ia2  *  ”2 

A2  2 
2  + 

which  is  (4.2).  iMs  proves  Theorem  4.1. 

More  qeneral  results  can  be  obtained  by  a  method  similar  to  that  of 
Widlund  (1978). 

Theorem  4. 2 

For  the  method  (2.16),  for  k  even  or  k  *  1, 
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(4.5) 


|rk|/|r°|  4  2 Pk/(1  +  P2k) 

and  for  k  odd 

(4.6)  |rk|/|r°|  4  2pk/(1  -  p2k) 
and 

(4.7)  |rk!/!r°!  <  4pk/((1  +  p2)(1  +  p2(k-1)))  , 

where  p  =  A /(/  i  +  A2  +  1).  Note  that  (4.6)  is  a  better  estimate  than  (4.7) 
when  A  is  large  and  k  is  large,  (4.7)  is  better  when  A  or  k  is 
small.  Note  that  (4.5)  for  k  =  1  and  k  =  2  gives  the  same  result  as 
Theorem  4. 1 . 

Proof  of  Theorem  4.2 


Bv  (3.4)  and  (3.7) 


where  z  is  in  the  span  of 
z  =  Pk(I+S)r°  where  Pk(A) 
Thus 


.  k  k.  .  k  0. 

(r  ,  r  )  *  (r  ,  z  +  r  ) 

(I+S)r°, . . . (I+S)rk-1 .  Now  by 


is  a  polynomial  of  degree  k 


(2.16) 
with  Pk(0) 


0. 


|rk|  =  (rk,  0k(I+S)r°) 

where  ( A )  is  a  polynomial  of  degree  k  and  0^(0)  =  1.  By  the  spectral 
mapping  theorem 


|rk|  4  min  max  ( 0.  ( 1  +  U)  I  |r°| 
0k(0)=1  P60(S) 


where  o(S)  is  the  spectrum  of  S.  Since  S  is  skew-symmetric  o(S)  is 
contained  in  the  imaginary  axis  with  -1  4  U/iA  4  1,  The  minimum  is  taken 
over  all  polynomials  0^ ( A )  of  degree  k  with  C^(0)  *  1.  As  does  Widlund 
(1*176),  we  take  the  particular  polynomial 


0k(1  +  U) 


Tk  ( M/i  A ) 
t  '(-i/i”A‘) 
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where 


T^(A)  =  cosh(k  cosh-1  A).  We  have  |Tk(U/iA)j  = 
|cos(k  cos  Vu/iA))!  <  1  and 


T^f-l/iA) 


-1 


=  cosh(k  cosh  (-1/iA))  =  cosh(k  loq 
=  j  ( (/l  +  A2  +  1  )k  +  (/l  +  A 


(Ay  +  A2  +  1 )  A  1 ) 


2  1)k)A-k 


1  2k  -k 

j  (1  +  P  )P 


for  k  odd  and 


t Tk ( — 1 /iA ) |  =  | sinh (k  loq{ 


/iT 


Az  +  l)A-1)| 


*  1  ((/l  +  A2  +  1)k  -  (A  +  A2 


1)k)A-* 


-1(1  -P2k)p"k 


for  k  even. 

This  proves  (4.5)  for  k  t  1  and  (4.6).  Inequality  (4.5)  for  k  =  1  is 

just  (4.1),  and  (4.7)  follows  from  (4.5)  and  (4.1).  This  proves  Theorem  4.2. 

k  v 

The  above  theorems  qive  results  on  the  error  vectors  e  since  r  = 


(I+S)eK  and  so 


(4.8) 


|ek|  <  | rk |  <  Ay  +  A2  |ek|. 
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5.  The  case  with  general  symmetric  part 

In  this  section  we  discuss  the  more  general  case  when  P,  the  symmetric 
part  of  A,  is  not  the  identity.  In  this  case  the  system 
(5.1)  Ay  =  (P  +  Q)y  =  c 

can  be  transformed  to  one  of  the  form  (1.1)  by  setting 


(5.2)  x  =  P1/zy  ,  b  =  P-1/2c  ,  S  =  p"1/2qP“1/2  . 

The  algorithm  (2.16)  can  then  be  used  on  the  resulting  system.  It  is,  of 
course,  more  convenient  to  work  with  the  original  matrices  P  and  Q  than 
with  S,  thus  we  consider  the  system  (2.16)  in  the  original  variables.  We 
have 


a) 

k+1 

w 

k  k 

*  w  +  a,  p 
k 

(5.3) 

b) 

k+1 

r 

k  k 

=  r  -Vp 

c) 

k+1 

P 

-1  k+1  0  k 

-  P  r  -  V 

where 

(5.3) 

d) 

“k  = 

(r  ,  P  r  )/  ( Ap  ,  p  Ap  ) 

and  8.  =  1  -  a.  . 
k  k 

The  vector  w° 

is  arbitrary. 

-  0  „-1  0  0 
and  p  =  P  r  where  r  =  c 

vectors  w*  converge  to  y  the  solution  of  (5.1).  Because  of  the  necessity 
-1  k+1  -1  k 

of  computing  P  r  and  P  Ap  the  method  is  applicable  in  practice  only 
when  the  solution  of 

Pz  =  d 


can  be  obtained  easily.  This  restriction  applies  also  to  the  method  of  Concus 

—  1  1(4. 1 

and  Golub  (1976)  and  Widlund  (1978).  Note  that  P  r  can  be  obtained  by 


-1  k+1  -1  k  -Ik 

Pr  *P  r  -a^p  Ap 


so  only  one  inversion  of  P  is  required  per  iteration  step. 
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The  orthogonality  results  for  the  algorithm  (5.3)  are  easily  deduced  from 


the  results  for  (2.16).  We  state  these  results  for  completeness. 
Theorem  5. 1 

For  the  algorithm  (5.3)  the  following  relations  hold. 

(A p* ,  P  ^p-’)  *  0  for  i  +  j 
(P-1rk,  AP“1r‘')  =  0  for  j  <  k 


,  -1  k 
(F  r  , 

Ap1)  s  0 

for  j  <  k 

.  k+1 
(P 

K 

N 

O 

for  k  >  0 

,  k  -1  j  .  k  -1  0  .  >  . 

(r  ,  P  rJ  )  =  (r  ,  P  r  )  for  j  <  k  . 

Theorem  5.2 

For  the  algorithm  (5.3)  the  following  estimates  holds. 


kk|  .,/k0!  <  2pk/ (  i  +  p2k) 

p  p 

for  k  even  or  k  *  i 

|rk|  _i  <  2pk/(1  "  p2k) 

P  P 

for  k  odd 

and 

| rk |  /|r°|  <  4pk/( ( 1  +  P2)(1  +  P2(k"1>)) 

p“  p~ 

for  k  odd 

where  |rk|2  =  (r\  P_1rk),  P  *  A/(/ 1  +  A2  +  1)  and  A  is  the  spectral 

P~' 

radius  of  P *1/2C)P’1/2. 

The  proofs  of  these  results  follow  from  Theorems  3.1,  3.2,  and  4.2. 
Corresponding  to  (4.8)  we  have 

(5.4)  (ek,  Pek)  4  (rk,  P_,rk)  <  +  A2  (ek,  Pek) 

which  in  conjunction  with  Theorem  5.2  gives  estimates  for  the  error. 
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6.  Numerical  Experiments 


The  performance  of  the  qeneral  alqorithm  (5.3)  depends  significantly  on 
the  means  of  inverting  the  positive  definite  matrix  P.  Therefore,  to  test 
the  algorithm  it  seemed  best  to  study  only  the  basic  algorithm  (2.16)  i.e. 
with  p  being  the  identity.  A  FORTRAN  computer  code  was  written  which 
implemented  the  algorithm  (2.16).  The  program  was  run  on  the  UNIVAC  1100 
computer  at  the  Madison  Academic  Computing  Center  using  single  precision 
arithmetic  which  has  about  seven  digits  of  accuracy. 

The  computer  code  also  implemented  the  alqorithm  of  Concus  and  Golub 
(1976)  and  Widlund  (1978).  We  will  refer  to  this  alqorithm  as  the  CGW 
algorithm.  The  CGW  algorithm  and  (2.16)  were  run  simultaneously  but 
independently. 

The  matrices  S  used  in  the  experiments  were  banded  skew-symmetric 

matrices  whose  non-zero  elements  were  generated  randomly.  For  0  <  i  -j  <  m, 

was  a  randomly  generated  floating  point  number  in  the  interval  [-<$,<5] 

for  some  number  6,  0  <  6  <  1,  also  =  -Sij'  otherwise  was  zero. 

For  all  the  numerical  experiments  x^1  =  0,  with  the  exact  solution  being 

Jc  0  _5 

=  1.  The  algorithm  (2.16)  was  stopped  when  |r  |/|r  |  <10  . 

Table  1  displays  the  results  from  several  experiments.  N  is  the  number 
of  unknowns,  m  gives  the  size  of  the  band  and  6  is  the  range  of  random 
numbers  for  S  as  described  above.  A  is  the  spectral  radius  for  S  as 
computed  by  EISPACK  routines  (Smith  et  al.  1976).  The  fifth  column  gives 
I,  the  number  of  iterations  required  for  convergence,  and  the  sixth  and 
seventh  columns  give  the  values  of  | e1 J / 1 e° |  for  both  the  alqorithm  (2.16) 
and  the  CGW  algorithm.  The  eighth  column  contains  the  quantity  appearing  on 
the  right-hand  side  of  the  estimate  (4.5),  for  k  *  I. 
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It  was  found  that  the  two  algorithms  converged  at  about  the  same  rate 
with  (2.16)  having  slightly  smaller  values  for  the  norms  of  the  error  and 
residual  vectors.  The  similarity  is  not  surprising  since  the  error  estimate 
(4.5),  usinq  (4.8),  is  similar  to  that  qiven  by  Widlund  (1978)  for  the  CGW 
alqorithm.  For  the  algorithm  (2.16)  it  was  found  in  all  the  experiments  that 
the  norm  of  the  error  decreased  monotonies lly.  This  was  different  than  the 
CGW  algorithm  for  which  the  norm  of  the  error  usually  did  not  decrease 
monotonically  for  the  first  several  iterates.  However,  for  the  CGW  algorithm 
the  even  and  odd  iterates  do  give  monotonically  decreasing  values  for  the  norm 
of  the  error,  widlund  (1978). 

The  estimate  (4.5)  is  seen  to  give  a  good  approximation  of  the  behavior 
of  the  alqorithm  and  can  be  used  to  give  a  good  estimate  of  the  number  of 


iterations  required  for  convergence 
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